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It is shown how to explicitly coarse-grain the microscopic dynamics of the rule-based Vicsek 
model for self-propelled agents. The hydrodynamic equations are derived by means of an Enskog- 
type kinetic theory. Expressions for all transport coefficients at large particle speed are given. The 
phase transition from a disordered to a flocking state is studied numerically and analytically. 



Pattern formation and collective motion in systems of 
self-propelled objects are fascinating phenomena which 
have attracted much attention. Systems of interest in- 
clude animal flocks [l| , chemically powered nanorods Q , 
and actin networks driven by molecular motors [3| ■ Theo- 
retical studies of these systems are usually based on phe- 
nomenological transport equations. In most cases, the 
equations are postulated by means of symmetry argu- 
ments, which define only the general form of the terms 
but leave their coefficients undetermined. 

One goal of this Letter is to provide a system- 
atic derivation of all relevant coefficients for the two- 
dimensional Vicsek model (VM) of self-propelled parti- 
cles In the VM, pointlike particles are driven with 
constant speed. At each time step, a given particle as- 
sumes the average direction of motion of its neighboring 
particles, with some added noise. As the noise amplitude 
decreases, the system undergoes a phase transition from 
a disordered state, in which the particles have no pref- 
ered global direction, to an ordered state, in which the 
particles move collectively in the same direction. This 
long-range order motivated renormalization group stud- 
ies by Toner and Tu ;5|. They found that the stabiliza- 
tion of the ordered phase is due to the nonzero speed 
of the particles, allowing two originally distant particles 
to interact with each other at a later time. The phase 
transition was originally thought to be continuous [J] but 
recent numerical work [fj indicates that the transition is 
discontinuous with strong finite size effects. There are 
few analytical studies on this transition 0, [|| . They do 
not treat the original VM but simple models related to 
it. For example, Bertin et al. [7[, study a model with 
simplified interactions and a continuous time dynamics 
by means of a Boltzmann equation. 

Numerical simulations of the VM 0, @] show localized 
high-density structures, for which a Boltzmann descrip- 
tion, which is restricted to low densities, is not sufficient. 
Enskog's proposal to generalize the Boltzmann equation 
to dense gases was a major milestone in kinetic theory. In 
this Letter, it is shown how an Enskog-type equation with 
genuine multi-body collisions can be obtained for the VM 
and how this can be used to rigorously derive hydrody- 
namic equations. In addition to the terms postulated by 
Toner and Tu the derived equations contain several 
new relevant terms which describe an intricate coupling 



between density and order parameter gradients. The co- 
efficients of all terms, compatible with the symmetries 
of the system, are calculated explicitly in third order of 
a gradient expansion. The new kinetic equation is used 
to determine the mean-field phase diagram of the VM, 
which agrees well with direct numerical simulations but 
disagrees with the results of a related continuous time 
model 0] ■ This shows the importance of explicitely tak- 
ing the discrete time, rule-based nature of the VM into 
account. The derived hydrodynamic equations are ap- 
plied to study the stability of a homogeneous flocking 
state against spatio-temporal perturbations. I discuss 
how an instability at the onset of collective motion can 
change the appearance of the phase transition from sec- 
ond to first order. Predictions for the system size where 
this change is expected to happen, are given. 

In the VM, a system of N pointlike particles with 
continuous spatial coordinates rj(i) and velocities v,(t) 
evolves via two steps: streaming and collision. During 
a time step r, particles stream ballistically: x^(i + r) = 
Xj (t) + TVi (t) . The magnitude of the particle velocities is 
fixed to vo. Only the directions 6i of the velocity vectors 
are updated in the collision step: a circle of radius R is 
drawn around a given particle and the average direction 
6i of motion of the particles within the circle is deter- 
mined according to 0j = arctan[^"sin(0j)/53™ cos (@j)]- 
The new directions follow as 9i (t + r) = 6j (t) + , where 
£i is a random number chosen with uniform probabil- 
ity from the interval [—77/2,77/2]. Since explicitly coarse- 
graining the dynamics of the VM is difficult, in previ- 
ous work Q, I have first validated the formalism on a 
simpler equilibrium model .10] which shares essential fea- 
tures with the VM. The kinetic formalism starts with the 
Liouville equation for the N-particle probability density 
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where XW = ( Xl) x 2 , . . . , x N ), 0W = (0 1 ,d 2 , . . . ,0 N ), 
and S(x) = J2m=-oc + 27rm) is the periodically 
continued delta function. The velocities = 
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(vi, V2, vjv), are given in terms of angle variables, 
v i = vq(cos 9i, sinflj). The collision integral contains 
integrations over the pre-collisional angles 9j. Assum- 
ing that the particles are uncorrelated prior to the col- 
lisions, the probability distribution can be expressed as 
a product of identical one-particle probability distribu- 
tions: P(6>W,XW) = nii-PiO^Xi). Thi s approxi- 
mation of molecular chaos is valid at moderate and large 
noise strength rj and when the mean free path (mfp) is 
large compared to the radius of interaction R. Here, the 
mfp is defined as the distance a particle travels between 
collisions, tvq, and is density-independent due to the dis- 
crete nature of the dynamics. Multiplying Eq. (TTJ) by 
J2i — v i)<5( x — x i) an d integrating over all particle 
positions Xi and angles Qi , yields in the large iV-limit Q , 
a kinetic equation for the one-particle distribution func- 
tion, f{9,x,t) = NP 1 {9,x,t), 
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where M/j(x, t) = J R p(y,t) dy is the average number 
of particles in a circle of radius R centered around x. 
The local particle density p is given as a moment of the 
distribution function, p(x,t) = J W f(9,x,t) d9; = 
J R ... <ix2 c?X3...(ix n denotes the integration over all posi- 
tions, 71—1 particles can assume within the interaction 
circle; (■■■)§ = J * ...d6id02...dO n is the average over all 
pre-collisional angles of n particles in the interaction cir- 
cle. Since particles in the VM have zero volume, there 
is a non-zero probability that a large number of parti- 
cles can be found in the collision circle of a given parti- 
cle. This leads to the unusual structure of the collision 
integral in which every term in the sum accounts for a 
n-particle collision. For example, the n = 4 term involves 
the product of four distribution functions and describes a 
four body collision. Interactions between particles which 
are not at the same position but a distance < R apart 
are explicitely taken into account by Eq. (|2|). This leads 
to collisional momentum transfer which is a key feature 
of the Enskog equation and not included in Boltzmann- 
type equations. Hence, Eq. ©, can be interpreted as an 
Enskog-like equation for pointlike particles with discrete 
time evolution; it remains valid even at infinite density. 

Let us first consider a spatially homogeneous system 
and study stationary solutions of Eq. ([2]) . This amounts 
to solving the fixed-point equation fo(9) = C(/o) for the 
stationary distribution function /o, where C denotes the 
r.h.s. of Eq. ([2]). It can be easily checked that the con- 
stant distribution /o = po/2ir is a fixed-point at any noise 
and average density, po — N/A, where A is the area of 
the system. This solution corresponds to the disordered 
phase, where all velocity directions occur at equal prob- 




Numerics: Nagy et al. 
▲ Theory: Aldana et al. 
□ Numerics 



a) 



M 




FIG. 1. a) The critical noise rjc as a function of the average 
number of collision partners, Al — pouR 2 , and the prediction 
of Eq. (35) for large vo from Ref. u\ , (dashed line), in com- 
parison with results from Refs. [1, \ElLM- b) Real part of the 
growth rate, ujr, of a small longitudinal perturbation of the 
ordered state versus dimensionless wave number at M — 5, 
very close to the threshold, (go —rj) /gc = 0.00057. The insert 
shows a lower and an upper bound for the crossover length 
L* (in units of the mfp) beyond which the phase transition is 
expected to become discontinuous. 

ability. Below a critical noise rjc(po) there exists another 
fixed-point solution which breaks rotational symmetry. 
It has a maximum at some arbitrary angle 9 and de- 
scribes ordered motion into this direction. The critical 
noise follows from the condition A = 1, with 
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Here, Mr is equal to nR 2 po and 9 is the average angle de- 
fined above Eq. (fTJ. The fixed-point equation was solved 
numerically for r\ <r\c- The solution approaches a cosine 
with vanishing amplitude when r\ approaches the critical 
noise. By means of a Fourier cosine series in 8 — 9 the 
behavior at the critical point was extracted analytically. 
The order parameter, defined as the amplitude g\ of the 
first non-trivial Fourier coefficient, is found to behave as 
gi oc \Jr\c — V- Thus, the order-disorder transition ap- 
pears to be continuous with the mean-field critical expo- 
nent of 1/2. Fig. la) shows the calculated phase diagram 
(solid line). Evaluating Eq. ([3]) in the low density limit 
gives rjc oc R^fp. This scaling with the square root of the 
density agrees with previous numerical [J] and theoreti- 
cal results 0, 0|- However, there is no dependence of the 
critical noise on the particle speed in the large mfp limit, 
which is consistent with numerical simulations of the VM 
0, EH but disagrees with the scaling r\c oc V pRvo for 
p — >• of the continuous model of Ref. [7[ . The dashed 
line in Fig la) shows that thephase diagram of this model 
(obtained from Eq. (35) in [7| with Vqt/R = 5) does not 
describe the VM. Evaluating Eq. ([3]) in the infinite den- 
sity limit yields rjc — > 2tt. In order to see whether the 
homogeneous ordered state is stable under time evolu- 
tion, I derive the hydrodynamic equations by means of a 
Chapman- Enskog expansion 0, [l2| ■ The basic idea be- 
hind this expansion is to take the local stationary state 



3 



as a reference state and expand around it in powers of 
the hydrodynamic gradients. To systematically account 
for these gradients a dimensionless ordering parameter e 
is introduced, which is set to unity at the end of the cal- 
culation. The procedure starts with a Taylor expansion 
of the l.h.s of Eq. (J2|) around (0, x,i). The spatial gradi- 
ents that occur are scaled as d a — > ed a , and multiple time 
scales U are introduced in the temporal gradients. These 
time scales describe different physical processes, for ex- 
ample, in regular fluids, the time scale proportional to e 
describes convection. For the VM, this is expressed as 
d t = d to + ed tl + e 2 d t2 .... 

Expanding the distribution function and the collision 
integral in powers of e, / = /o + efi + f 2 f 2 + and 
C = C + eCi + e 2 C 2 + inserting into Eq. ©, and 
collecting terms of the same order in e leads to a hi- 
erarchy of evolution equations for the Due to the 
absence of momentum conservation and Galilean invari- 
ance this set of equations is dramatically different from 
the usual one. It is not a priori evident whether the 
scaling ansatz for the time derivatives is correct. How- 
ever, it turns out that this choice avoids any incon- 
sistencies if additionally the expansion of the distribu- 
tion function / is identified as an angular Fourier series 
with /o(x, t) = p(x, t)/2n and, for n > 0, / n (x, 9,t) = 
[a n (x, t) cos (n0) + &„(x, t) sin (n0)] /ttvq. 

Many moments of the collision integral such as 
(v x v y C 2 ) — J q 2jt v x v y C 2 d0 are required in the Chapman- 
Enskog expansion. For simplicity, these moments are 
evaluated in the limit of large mfp, tvq ^> R. This in- 
volves solving the following four integrals, 
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The average angle is 



is given by vPi 
cos sin cos 0i sin 02 , ^3 = 
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a function of the angles 0i , 02, . . . 6 n . 

We seek a hydrodynamic description of the first two 
moments of /, namely the particle density p = J Q n f d9 
and the macroscopic momentum density vector w — 
(w x ,w y ), w = J Q 27r vfd9. Inserting the Fourier represen- 
tation of / into these moments shows that the first order 
coefficients are given by the momentum density, a\ = w x 
and bi — w y . Multiplying the hierarchy of evolution 
equations by powers of the microscopic velocity vector 
v = (v x ,v y ) and integrating over gives a set of equa- 
tions for the time development of the density and the 
moments and bi. This analysis is performed in the 
vicinity of the critical point, |A — 1| <C 1, in order to sig- 
nificantly simplify the consistent closure of the hierarchy 
of moment equations, see 
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For simplicity, all equations are rescaled by expressing 
time in units of r and distances in units of the mfp, tvq, 
which also makes p and w dimensionless. After straight- 
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TABLE I. The transport coefficients hj, qj and kj, defined in 
Eq. (J6]), are expressed as functions of F, S, p, q, see Eq. (|8]). 

forward, but tedious, calculations one obtains the con- 
tinuity equation dtp + d a w a — 0, and a rotationally- 
invariant equation for the momentum density, 

dtvS + V H = -Wp+(\- Vjw + Qx -w + Q 2 -Vp (5) 

with b = (3 — A)/4. The momentum flux tensor H and 
the tensors Qi, Q2, 
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are given in terms of five symmetric traceless tensors 

£li, a p = d a wf) + dpw a - 8 af) d 1 w 1 
^2, a /3 = 2d a dpp - 8 a pd 2 p 
fl 3 , a /3 = 2w a wp - S a/3 w 2 

^4,a/3 = W a dpp + WfjdaP - 8 aP w 1 d 1 p 

^5, a/3 = 2{d a p){dpp) - S a0 (d 7 ) 2 . (7) 

The tensor f?i is the viscous stress tensor of a two- 
dimensional fluid. The transport coefficients in Eq. © 
are given in Table HI They depend on the following vari- 
ables, 
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where 7 is the ratio of the interaction radius to the mfp, 
7 = R/tvq. Eq. ((5} is consistent with the one postulated 
in Ref. [j| but contains additional gradient terms. It has 
a homogeneous flocking solution: w = wq n and p = pq. 
The amplitude of the flow is given by ivq = y(l — A) /q%. 
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In order to study the spontaneous onset of collective 
motion, a perturbation around this state is considered, 
p{x,t) = p + dpe lk - x+UJt , u?(x,t) = w n + (5we ik ' x + wt , 
and Eq. ([5]) is linearized in dp and Sw. The characteristic 
equation for the growth rate uj(k) describes three possible 
modes. I found that in a small window, 7/5 < r) < f]c, di- 
rectly below the onset of flocking, one of the longitudinal 
modes is always unstable against long wavelength pertur- 
bations: the real part of uj is positive for < k < ko as 
shown in Fig. lb). A similar instability was reported by 
Bertin et al. Q- Chate et al. Q found numerically that 
the order/disorder transition is discontinuous for system 
sizes L larger than the crossover length L*. Assuming 
that the long wave instability is the reason for this fi- 
nite size effect, I calculated the largest value of fco within 
the narrow instability window at constant density, k* , in 
order to obtain a lower bound for L*. Plotting 2ir/k* 
gives the lower curve in the insert of Fig. lb). An up- 
per bound was obtained by determining the wave number 
kmax where the growth rate has the largest value inside 
the instability window. The upper curve in the insert 
shows 2it/k max as a function of density. The minimum 
around M ss 2 and the divergences at small and large 
densities are consistent with numerical results @. 

To see what happens to a growing perturbation be- 
yond the linear instability, the continuity equation and 
Eq. © were integrated on a L X L lattice with peri- 
odic boundaries by means of a predictor-corrector scheme 
13]. These simulations confirmed that the ordered 
phase is stable for small system sizes L < 2-7r/fco- For 
slightly larger system sizes one observes a stable, inho- 
mogeneous steady state with a global order parameter, 
(w) = J w dx/L 2 , larger than the amplitude of the homo- 
geneous state, w . Finally, for much larger system sizes, 
it turns out that the system is both linearly and nonlin- 
early unstable for rjs < 77 < rjc- Longitudinal pertur- 
bations grow without bound; they do not lead to stable 
solitons as suggested in Ref. [7J. However, direct simula- 
tions of the VM at large mfp do show solitary structures 
such as traveling high-density bands in a window just be- 
low the transition jy, |Tl| . At lower noise these structures 
disappear. Identifying this "solitary" window with the 
instability window, its size can be predicted by the cur- 
rent theory which takes all the details of the VM such 
as multi-body interactions into account. However, inside 
this window, the hydrodynamic equations are driven out 
of the range of their validity and are not suited to describe 
solitons. Nagy et al [4| did not see high-density bands at 
small velocity vq. To treat this limit of small mfp theoret- 
ically, one has to abandon the molecular chaos approx- 
imation i.e. go beyond the mean-field approximation, 
which is outside the scope of this paper. 



In summary, a first-principle derivation of the hydro- 
dynamic equations of the VM by means of a novel ki- 
netic theory is presented and a stability analysis of the 
resulting equations, Eq. ([5]), is performed. The mean- 
field phase diagram for arbitrary density is calculated. It 
agrees within a few percent with simulation results and is 
shown to be independent of the particle speed in the large 
mfp limit. It is also shown that the continuous theory of 
7] fails to reproduce the phase diagram of the VM and 
that one has to explicitely incorporate the discrete time 
dynamics and genuine multi-body interactions in order 
to achieve agreement. The theory presented here is con- 
sistent with numerical studies [!, |(| , and suggests the fol- 
lowing picture of the nature of the flocking transition in 
the large mfp limit considered here: At 7/ = r\c a homo- 
geneous ordered state bifurcates continuously from the 
disordered state. At the threshold, this state is unstable 
to longitudinal, long wavelength fluctuations. Perturba- 
tions from a large range of wave numbers k < ko become 
unstable, already in close vicinity to the threshold. The 
transition appears to be continuous in small systems but 
becomes a discontinuos transition in large systems due 
to the emergence of density waves which abruptly in- 
crease the global order parameter. An estimate of the 
system size L*, above which the discontinuous nature of 
the transition is expected to emerge, is given. This length 
is found to diverge at small and large densities, consistent 
with numerical results. 
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